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Abstract 


R. Ewing, O. Liev, R. Lazarov and A. Naumovich in [1] proposed a finite volume 
discretization for one dimensional Biot poroelasticity system in multilayer domains. 
Their discretization and exact solution are invalid. We derive valid discretization 
and exact solution. Finally, our numerical solution is compared with known exact 
solution in discrete Lg norm. 
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1 Introduction 


The presence of a moving fluid in a porous medium affects its mechani- 
cal response. At the same time, the change in the mechanical state of the 
porous skeleton influences the behavior of the fluid inside the pores. These 
two coupled deformation-diffusion phenomena lie at the heart of the theory 
of poroelasticity. More precisely, the two key phenomena can be summarized 
as follows: 


1. fluid-to-solid coupling: occurs when a change in the fluid pressure or the 
fluid mass induces a deformation of the porous skeleton. 

2. solid-to-fluid coupling: occurs when modifications in the stress of the 
porous skeleton induce change in the fluid pressure or the fluid mass. 


In accordance with these two phenomena, the fluid-filled porous medium acts 
in a time-dependent manner. Indeed, suppose that the porous medium is com- 
pressed. This will result in an increment of the fluid pressure inside the pores 
and consequent fluid flow. The time dependence of the fluid pressure will in- 
duce a time dependence of the poroelastic stresses, which in turn will respond 
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back to the fluid pressure field. The earliest theories, which is related to Terza- 
ghi, accounted for the fluid-to-solid coupling only. In this case, the problem is 
mathematically much easier. This kind of theory can model successfully some 
of the poroelastic processes in the case of highly compressible fluids such as 
air. However, when one deals with slightly compressible or incompressible 
fluids, the solid-to-fluid coupling cannot be neglected since the changes in 
stress field can influence significantly the pore pressure. Maurice Biot was 
the first who, by means of phenomenological approach, developed a detailed 
mathematical theory of poroelasticity which successfully incorporated both 
basic phenomena mentioned above. In this paper, assumption of only ver- 
tical subsidence is invoked and this leads to the one dimensional model of 
poroelasticity. We consider a finite volume discretization for one dimensional 
Biot poroelasticity system in multilayer domains. For stability reasons, stag- 
gered grids are used. The discretization takes into account discontinuity of 
the coefficients across the interfaces between layers with different physical 
properties. 


2 Biot model in one dimension 


In one dimension, the domain of consideration (2 is an interval (0, Z) where 
the boundary I is {0, L}. The Biot model, which describes poroelastic process 
in (2 can be written as a system of partial differential equations for the 
unknown fluid pressure p(x, ¢) and displacement of the porous medium u(z, t) 
consisting of the equilibrium equation and the diffusion equation 


0 Ou Op _ 
— gy (A+ 2H) a) + a = 0, x © (0,2), t € (0,7), 
a) _ Ou O «Op, 
54 PPP | Fa? Be aoe x € (0,L), t€ (0,7), 


where A (dilation moduli) and yz (shear moduli) are Lame coefficients of the 
porous medium. Here ¢, 6, «, 7 and q(x,t) are porosity of porous medium, 
compressibility of the fluid, permeability of the porous medium, viscosity 
of the fluid and source term, respectively. We define stress tensor and fluid 
velocity, respectively by the following relationships 


Ou k Op 


eS eas 1 Ox" 


In classical formulation, the one-dimensional Biot model describes, fluid flow 
and skeleton deformation caused by the constant vertical load applied on the 
top of column of soil, which is bounded with rigid and impermeable bottom 
and lateral walls, and a top wall which is free to drain. The following boundary 
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and initial conditions supplement this model 


O 
p=0, O+2)— =-s, at r=0. 
Ox 
This means that the upper boundary is free to drain and a load with the 
value so is applied to it. Also 


means that the variation in water content is zero at the beginning of the pro- 
cess. Now, consider the case when the porous medium is not homogeneous but 
has a layered structure, each layer being characterized by different porosity, 
permeability and Lame coefficients. For the simplicity of presentation, let us 
restrict ourselves to the case of only two layers. In the case of the considered 
two-layered medium, coefficients of the governing equations are discontinuous 
across the interface € 


x >, M2 u>&, 


n(x) ={ LSE, o- {3 r<E, 


K2 x> €, 2 > €. 


Assuming a perfect contact, the interface conditions look as follows 
[u] = 0,  [p] = 0, (1) 


which express continuity of the displacement and of the fluid pressure across 
the interface. Also 


[Ss]}=0, [V]=0, (2) 


which means continuity of the stress of the porous skeleton and continuity of 
the fluid flux, respectively. In the formulae (1) and (2), we have 


[a] = |a=€+0 —q le=€—0; 


where q is a symbol for quantities u, p, S and V. As it is shown in [3], the set 
of interface conditions (1) and (2) can also be derived directly from the Biot 
equations if they are written for a general inhomogeneous medium. Now, the 
following dimensionless dependent and independent functions are introduced 
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ee” oo, PO og, ge 
L bo mol? 7 89" 7 
K 
A + 2 n Ln 
= .——__—,, = p=, = Ao +2 ; vt = #). 
OS og? me O bB(Ao + 2H0), f(x, t) aon q(x, t) 
10 


Then, the governing equations together with the boundary, initial and inter- 
face conditions can be transformed to dimensionless form 


Oo, Ou Op 
= 1 cg 

ag an) + Bp =O «FED, €€ (0,7), 
0 Ou 0, OD. _ 
7 (OP Fa) 5g Hq) = Fist), e401), Le (0,7), 
pot = 1, p=0, at x=0, te [0,7], (3) 

=0, oe a at x=1, te (0,7), 

Ox 
fps at t=0, x€ (0,1), 
Ox 

[u] = 0 we =0 [p] =0 i a t =€, te [0,7] 
U = + Ox ~~ ? Pp > ? Ox i ? a ed ? ? 


Further, the possible discontinuities of the dimensionless coefficients at « = € 
are distinguished 


(a) = {¥ eat ne) ={ oe, 


V2 &> €, K2 a> €, 


a(t) = {8 wee, 


dg «£>€. 


For the convenience of the theoretical analysis, the problem (3) is transformed 


into a problem with homogeneous boundary conditions, by the following sub- 
stitution 


L il 
u(a,t) = u(#,t) -— —a4+-. 
Vv v 


According to this substitution, problem (3) is reformulated as follows: 
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QO, Ou Op _ 

aa’ Da? i Ox 0, ee (0,1), Pe (0,7), 
0 _ Ou ip OD. 
74 oP + Fa? 5g an) ~F(e 4), “= (0,1), £€ (0,7), 
Ou 
vx = 0, p=0, at x=0, te [0,7], (4) 
u=0 noe =o, at x=1, te (0,7), 

Ou il 

ap+ a =o at t=0, x €(0,1), 
ful=0, WS) =0, pl=0, [2)=0, at 2g, t€ (0,7) 
Ul = ’ YA mal o] D ~~ ’ Y Bas ~~ ? a ed ;] ’ ‘ 


Due to the complexity of the Biot system, analytical solutions in closed form 
are available only in very special cases. Certainly, the situation gets compli- 
cated in the case of inhomogeneous porous media. The choice of the numerical 
method for the discretization of the poroelasticity system is not obvious. The 
finite element method currently dominates in solving poroelasticity system, 
especially when dealing with complex domains (see [4] for further details). 
Although finite element methods can be applied to the interface problems, 
however, they usually work on grids which resolve the interfaces. Hence this 
fact leads to that imposes certain restriction on the method. Moreover, even 
when the grids resolve the interfaces, standard finite element methods do not 
provide good approximation for the flux variables. On the other hand, there 
is variety of successful finite difference and finite volume approaches, where 
the interfaces are allowed to cross the grid cells (see [5]). 


2.1 Grids and notations 


For the interval (0,1) and N > 1, we define stepsize h in the following form 


h:= z : 
2N—1 


To overcome stability difficulties, which often arise when the discretization of 
the Biot model is done on the collocate grids, the use of staggered grids was 
proposed in [2], [7]. Two different spatial grids, @, to discretize the pressure 
equation and &, to discretize the displacement equation, are employed 

@p ={t,: 2; =th, i=0,1,..,N—1}, 

Oy ={Xi-0.5: i-0.5 = 2i—0.5h, i =1,2,..., N}. 


Further, the grids w, and w, are also used 
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Wp ={x; Ep, t= 1,2,...,N —1}, 
Wy ={®i-05 € Oy, +=1,2,...,N—1}. 

A grid in time with a stepsize 7 is also defined 
wr =f{tn: th=nt, n=1,2,...,M}. 


These grids as designed to represent the values of the pressure p at the grid 
points x; € W, and the values of the displacement u at the midpoints x;_9.5 € 
@, of the subintervals (x;_1,2;). According to these grids, position of the 
interface € could be represented in the form 


€ = %,,,,-0.5 + Oh, 


where 0 < ting < N is an integer and 0 < 6 < 1. Now, the following notations 
for discrete functions, defined on &, xX wr and W, X wr, respectively, are 
introduced 


U:=u" := uP := u(ai_o.5, tn), 
pi =p" = pi = p(xi,tn), 
p? :=op"*! + (1—a)p”, 


Moreover we use some notations for the first order forward and backward 
finite differences on a uniform mesh in the following form 


Dip, b) — play, t 
Px ee me a PA - ) 


vi,t) — p(a_i,t 
=e ist) vt — ) 


In a similar way we define 


3 t)— ~~ t 
Ug := Ue, = u(xito.5,t) - u(@i_o.5; ) 
Ug = Uzi = u(zi-o.s,t) . u(ti-1s,t) 


Finally, the finite differences in time are defined 


yer _ ue 
a See — si a : 
Uz = UL = Uy (Zs_0.5,tn) = = » 0.5 © Wu, 
n+1 n 
ba ee a _ Dd; Pi 
Pt = Pp = pi(&i, tn) = > Li EWy. 
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2.2 Finite volume discretization 


In order to approximate the differential problem (4) by finite volume method. 
Firstly the Biot equations are rewritten in the following way 


J i x€(0,1), te (0,7), 
Ox Ox 
6) O OV 2 
U 
5p oP i ag) + Ox = f(z,t), LE (0, 1), te (0, T]. 


Now, the first equation in (5) is integrated over the interval (a;_1, x;) 


PS ye OR es 
[ Td [- no (6) 


4-1 


and the second equation over the interval (20.5, Yi+0.5) 


Li+0.5 fa) Ou Zi+0.5 OV Zi+0.5 
/ 5p oP | 7a | ate | f(a, t)dz. (7) 


i—0.5 i—0.5 i—0.5 


Hence, in accordance with the interface conditions (1) and (2), some integrals 
from (6) and (7) can be rewritten as 


/ an ot = 9(x;) — S(x;_-1), / ——dz = V (xi+0.5) — V(xi-0.5), 
Vi-1 x Li-0.5 Ox 


Ly fe) Ti+0.5 0 
/ oF da = p(x) — plait); / SA dx = u(ai40.5) — u(ai-o.5)- 


Ox i-0.5 - 


Using the rectangular quadratic formula, we can write 


Li+0.5 fa) 0 Li+0.5 n+1 amt 
/ 5, (ap)da © in) | a(a)de me a;et—? 
x x 


LS a 
i-0.5 ot ot i-0.5 7 
where 
Li+0.5 
aj -| az dx (9) 
Li—-0.5 


In order to approximate the fluxes S(a) and V(x) in (8) in the grid points, 
with integrating the equation 


over the interval (x;_0.5,2ij+0.5) and the equation 
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over the interval (x;-1,2;) we will have the following integral equations 


Li+0.5 Li+0.5 xy Ly 
/ NOD = / aor / ve) dx = / op dx. 
Li-0.5 iy Li-0.5 Ox a1 © Li-1 Ox 


i- ir 


Then, by applying approximate formulae for integrals, we can transform these 
equations into the following form 


Lit+0.5 dx vy dx 
i Ty ak SS Ui —Ui-0.5) i-0. —~ dx = —(pi—Pi-1). 
sas) f way Uito.5—Ui-o.s, V(x-0 yf se) x (pi—pi-1) 


From these two formulae, approximate expressions for fluxes can be found 


Ui — U;_ asaya 
S(ai) ~ Si a , a a V (xi-0.5) ~ V; = —K oe, 


_ 1 ftte8 de yf de 
maf ay 8G [ @ | (10) 


After the substitution of approximate expressions for all the integrals into 
equations (6) and (7), weighted discretization in time with the weight pa- 
rameter o is applied. This procedure produces a finite difference scheme, 
which is a discrete analogue of the problem (4). The obtained finite differ- 
ence scheme is theoretically investigated and detailed convergence analysis 
is presented in [6]. Using non-index notations, this scheme for the discrete 
approximate solution u = ul!’ at point (a;~0.5,tn) € wy X wr and p = pi at 
grid point (x;,tn) € Wp X wr can be written as in the following form 


where 


Vy Z<E, 
V2 z>€. 


—fub+p$=0, c=295 G=1), tewr, (a) ={ 


— 1 (us), +ph=0, a <€, @=2,3,..,N—1), teuwr, 


S t HL ye + pS = 0, t1 <£, £5 >, (¢=2, 3, .,¥ — 1), tewr, 


— (us), +p =0, r_1>€, (¢=2,3,..,N—1), teur, 


(ap+uz)e— Poe + Soe =f", ios <& titos > &, @ = 1,2,..,N — 
2), tewr, 


(ap + Ux )t — K2(pZ) x a Li—0.5 > E, (i= 1,2,...,N — 2), tewr, 


(ap + Ue)e + Fp? = f%, © = an-1, (@ = N-1), t Cw, k(x) = 
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Ki ©N-1.5 S&, 
Ko EN-1.5 > €. 


po =90, un =0, teu, 
ap+U, = 4, f=2;€0,, (= 1,2,..,N =—1), t=, 


where coefficients a, « and v are calculated according to the formula (9), 
(10) and the right hand side f is defined as 


ee - / ee aa 


i—0.5 


3 Numerical results 


In this section, results of the numerical experiment are presented. Conver- 
gence of all unknowns of the system, i.e., w and p produced by the scheme 
(11) with respect to the exact solution of the continuous problem are shown. 
The numerical solution is compared to the known exact solution in discrete 
Ly norm, which is calculated according to the following form 


lewllas =h x. |titeee (ay Pace) = Teen Core ere | 
Li EG Wy 


where Wert and Wapp Stand for the exact and numerical solutions, respectively 
and w = {u,p}. In the following experiment, weight parameter is o = 0.5. 


Example 3.1 Suppose the following values of the parameters are used: 


87 
1 1 
Ay = 4, ? 
' 87 tan( 25) tan(222) 
a,=0, a2=0, f(a,t) =0. 


The position of the interface is at € = z. There is no exact solution of problem 
(4). Now, consider the following initial condition 


ipa <6 at t=0. 
Ox 


If we substitute the above condition in problem (4), the exact solution is as 
the following forms 
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cos(48") sin(#) exp(—0.25t) x <6, 
p(x, t) = 
sin(4) cos(4(1 — x)) exp(—0.25t) a> &, 
2cos(42*) cos($) exp(—0.25t) x < &, 
u(a,t) = 


cos _ . 
ae sin(47(1 — x)) exp(—0.25t) 


an 


x> €. 


3 


Note that the mesh size h is decreased in a way, preserving a constant value 
for the parameter @ in the expression € = 2x;~0.5 + 0h. The convergence 
results are given for times ¢ = 0.1 and t = 1. All numerical results are shown 
in Tables 1 and 2. In Figures 3.1(a-e) and 3.2(f-j), we have the convergence 


Table 1: Convergence in discrete Lz norm at the time ¢ = 0.1. 


RET Wlewlh llepll 
r 0.010978550724151 _[2.584356430578406 x 10-° 


£ 4.299657739862227 x 10—7 |1.565023743780048 x 10~© 
ign 5.255105309209073 x 10~® ]4.642021997862129 x 10-* 
= |5-263842780036844 x 10—19]2.411707866280830 x 10-7 
=e 3.408942895693747 x 10—11]6.152003882582387 x 10—§ 


Table 2: Convergence in discrete Lz norm at the time t = 1. 


- llewlh lleell 
6.010657532287119 x 10~® |2.069371443954429 x 10-> 


2.741582817269723 x 10~7 


9.979031969906696 x 10~7 


3.350803084864384 x 107% 


2.959883906284094 x 10~* 


3.356374342242523 x 10719 


1.537772829030507 x 10~7 


2.173637957624389 x 10~ 17 


3.922690864486794 x 10-8 


of displacement for given stepsizes h=0.1, 1/40, 1/160, 1/640 and 1/2560 at 
t=0.1 and t=1, respectively. Also, Figures 3.3(k-o) and 3.4(p-t) are prepared 
for representation of convergence of pressure for given stepsizes h=0.1, 1/40, 
1/160, 1/640 and 1/2560 at t=0.1 and t=1, respectively. In all of these figures, 
exact and approximate solutions are represented by continuous lines and 
broken lines, respectively. 
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Figure 3.1(a) 
i 


a 


h=0.1, t=0.1 


———— 
| | | | | i | 1 | 
O41 0.2 03 04 0.5 0.6 07 08 09 
x 
Figure 3.1(b) 
i 
h=1/40, t=0.1 
eg 
if ik 1 1 bt if 1 1 i 
O14 0.2 03 04 0.5 06 O7 08 og 


Figure 3.1(c) 
T 


h=1/160, t=0.1 


0.1 


0.2 


0.3 


0.4 


07 
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Figure 3.1(d) 
T 


h=1/640, t=0.1 


0.1 0.2 0.3 04 05 0.6 07 08 09 1 


Figure 3.1(e) 
i 


h=1/2560, 


0.1 0.2 0.3 0.4 05 0.6 07 08 09 1 


Figure 3.2(f) 


0.4 0.2 0.3 0.4 05 0.6 07 08 09 1 
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Figure 3.2(g) 
T 


i T i 
as = h=1/40, =1] 
Mise 
id if i L if i i L - 
01 0.2 0.3 04 05 0.6 07 0.8 0.9 2 
x 
Figure 3.2(h) 
T T T 
h=1/160, =1] 
ng 1 1 L Nd if i 1 i 
01 0.2 0.3 04 05 0.6 07 0.8 0.9 1 
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Figure 3.2(i) 


h=1/640, t=1 


0.1 0.2 0.3 0.4 0s 0.6 07 08 09 1 


67 


68 M. Namjoo, H. Atighi Lorestani 


Figure 3.2(j) 


h=1/2560, t=1 


0.41 0.2 0.3 0.4 05 0.6 0.7 08 09 1 


Figure 3.3(k) 
i 


rr. 


O41 02 03 04 0.5 06 O7 08 og 1 
x 
Figure 3.3(l) 
uf 
. h=1/40, t=0.1 


Finite volume method for one dimensional biot poroelasticity system... 


Figure 3.3(m) 
i 


0.1 0.2 0.3 04 05 0.6 07 0.8 09 1 


Figure 3.3(n) 
T 


™ h=1/640, t=0.1} 9] 
% 
es es J 
1 Nt i 1 L v4 1 1 i 
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Figure 3.3(0) 
T 


h=1/2560, t=0.1 


0.1 0.2 0.3 04 05 0.6 07 0.8 09 1 
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Figure 3.4(p) 
i 


neo, 1] 4 


0.1 02 0.3 04 0.5 0.6 07 08 09 1 


Figure 3.4(q) 
i 


= 
f™ h=1/40, t=1 7 


Nee 7 
1 1 1 1 1 1 1 1 1 
O41 0.2 03 0.4 Os 06 07 o8 0.9 1 
x 
Figure 3.4(r) 
i 
Pee: 
a h=1/160, t=1 
be 4 
1 1 1 1 1 1 1 1 1 
O41 0.2 03 0.4 Os 06 07 08 0.9 1 
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Figure 3.4(s) 
Mi 
7 h=1/640, t=1 x 
/ 4 
a 
if 1 i i L f i 1 L 
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x 
Figure 3.4(t) 
T 
~ h=1/2560, t=1 =] 
Xe 
see 7 
i L i i i i i i i 


0.4 0.2 0.3 04 0s 0.6 07 08 09 1 


4 Conclusion 


In this paper, we derived a finite volume discretization for the Biot system 
with continuous coefficients. In example 3.1, numerical solution was compared 
to the known exact solution in descrete Lj norm. In Tables 1 and 2 and de- 
rived figures, when the stepsize h becames smaller, we derived better results. 
In fact, our numerical experiments confirmed the theoretical considerations. 
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